Direct Numerical Simulation of Turbulent Heat Transfer Modulation in 

Micro- Dispersed Channel Flow * 
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Abstract The object of this paper is to study the influence of dispersed micrometer size particles 
on turbulent heat transfer mechanisms in wall-bounded flows. The strategic target of the current 
research is to set up a methodology to size and design new-concept heat transfer fluids with properties 
given by those of the base fluid modulated by the presence of dynamically-interacting, suitably- 
chosen, discrete micro- and nano- particles. 

We run Direct Numerical Simulation (DNS) for hydrodynamically fully-developed, thermally- 
developing turbulent channel flow at shear Reynolds number Ret = 150 and Prandtl number Pr = 3, 
and we tracked two large swarms of particles, characterized by different inertia and thermal inertia. 
Preliminary results on velocity and temperature statistics for both phases show that, with respect 
to single-phase flow, heat transfer fluxes at the walls increase by roughly 2% when the flow is laden 
with the smaller particles, which exhibit a rather persistent stability against non-homogeneous 
distribution and near-wall concentration. An opposite trend (slight heat transfer flux decrease) is 
observed when the larger particles are dispersed into the flow. These results are consistent with 
previous experimental flndings and are discussed in the frame of the current research activities in 
the field. Future developments are also outlined. 

PACS numbers: 



INTRODUCTION 



The problem of developing efficient heat transfer tech- 
niques for technological appHcations has become more 
and more important over the last decades due to the in- 
creasing demand of cooling in high heat flux equipments 
and to the unprecedented pace of component miniatur- 
ization |1|2|3| . Consider, for instance, teraflop computers 
and other electronic equipments Hke optical fibers, high 
energy density lasers and high power x-rays. These de- 
vices are required to operate with high precision and, 
at the same time, with minimum size. Such require- 
ments impose a challenge in terms of both device de- 
sign and thermal management, not only in the case of 
micro-scale applications but also for large-scale applica- 
tions such as transport vehicle engines, fuel cells and con- 
trolled bio-reactors [3]. Common air-based cooHng sys- 
tems have proven inadequate in high heat flux applica- 
tions and more efficient techniques for heat transferabil- 
ity are thus required. In particular, a quest for a fluid 
with high heat transfer capacity, characterized by the 
possibility of tunable thermal properties and also asso- 
ciated to low management/safety problems has started. 
One possibility is to use nanofluids, namely dilute Hquid 
suspensions of nanoparticulate solids including particles, 
nanoflbers and nanotubes, which are supposed to change 
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the heat transfer capabilities of the base fluid up to the 
target, desired amount. Nanofluids were flrst brought 
into attention approximately a decade ago, when their 
enhanced thermal behavior was observed with respect to 
conventional single-phase fluids such as water, engine oil, 
and ethylene glycol [3] . Speciflcally, due to the high ther- 
mal conductivity of metals in solid form, fluids containing 
suspended metal particles display signiflcantly enhanced 
thermal conductivity and speciflc heat capacity |4|5) com- 
pared to the conventional heat transfer liquids. For ex- 
ample, the thermal conductivity of gold at room temper- 
ature is more than 500 times larger than that of water 
and more than 2000 times larger than that of engine oil. 

The idea of increasing the effective thermal conductiv- 
ity of fluids with suspensions of soHd particles is not so 
recent since the flrst theoretical formulation of nanofluids 
as a new concept of heat transfer fluids was put forth by 
Maxwell [6] more than one century ago, then followed by 
several other experimental and theoretical studies, such 
as those of Hamilton and Grosser [7] and Wasp et al. 
[H]. As of today, however, a clearcut understanding of 
the modiflcations of the heat transfer mechanisms occur- 
ring in nanofluids is still to be produced. This is indi- 
cated by the failure of classical models of suspensions 
and slurries in predicting nanofluids behavior |3I9| ac- 
companied by the lack of alternative modelling strate- 
gies. Such lack may be possibly due to the number of 
investigations on the physical mechanisms which govern 
the heat transfer processes, still relatively low if com- 
pared to the corrispondingly large number of theoretical 
and experimental studies devoted to the problem of en- 
hancing heat transferability with nanofluids (see |3|9) for 
a review). A possible way to improve the understand- 
ing of the physical transfer mechanisms is to use accu- 
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rate and reliable numerical tools such as Direct Numeri- 
cal Simulation (DNS) and Lagrangian Particle Tracking 
(LPT), which may complement complex and costly ex- 
periments. DNS-based Eulerian-Lagrangian studies have 
been widely used for investigating mass, momentum and 
heat transfer mechanisms in turbulent boundary layers. 
In particular, previous DNS studies on turbulent particle 
dispersion in wall-bounded flows (see |10m2j among oth- 
ers) have proven their capability of predicting particle- 
turbulence interactions: these studies made it possible to 
perform phenomenological and quantitative analyses on 
the dispersion processes |13|14j . highlighting clustering 
and segregation phenomena. In the process of transfer- 
ring this research methodology to nanofluids, one ques- 
tion requiring clarification for clearcut expectations is to 
what extent the same approach is proper to solve the 
different new physics challenges imposed by the different 
phenomena and range of parameters. 

In this work, we propose to start approaching nanofiu- 
ids from a simplified simulation setting where micropar- 
ticles, rather than nanoparticles, are dispersed into the 
fiow. This simplification allows to single out complica- 
tions arising when the particles are very small and com- 
pare with the molecular scales of the fluid: inter-particle 
forces, wall effects, van der Waals forces, Brownian dif- 
fusion, etc. DNS has been used already to investigate 
on turbulent heat transfer in wall-bounded flows |15ll2T| 
and available studies provide systematic analyses of the 
Reynolds and Prandtl number effects on the heat transfer 
process. However, all these studies consider either single- 
phase turbulent flows |16|18ll2T| or turbulent flows laden 
with coarse particles |15|17j . A comprehensive analy- 
sis accounting for mass, momentum and heat transfer 
mechanisms all together and tailored for the specific case 
of micro- or nano-dispersed fluids is currently unavail- 
able as far as our knowledge goes. This is due to the 
non-trivial modehng issues, which of course reflect upon 
the complex interactions between the two phases. To 
elaborate, studying heat transfer modiflcations requires 
modeling particles as active heat transfer agents which 
interact both with the temperature field and the velocity 
field. Necessary energy and momentum coupling terms 
must be incorporated in the governing equations of both 
phases. This paper represents an effort toward a sys- 
tematic phenomenological study of turbulent heat trans- 
fer mechanisms in micro- and/or nano-dispersed fluids. 
Since the modeling of these fluids represents a largely 
unexplored field of research, this study involves substan- 
tial challenges due to the rich complexity of the involved 
physics. The main focus of the present paper is to ex- 
amine the modifications produced by solid inertial par- 
ticles on the temperature fields of both fiuid and parti- 
cles. First, the numerical methodology that we use to 
investigate on the problem will be described; then pre- 
liminary statistical results will be shown and discussed 
in the Hmit of hydrodynamically fully-developed, ther- 
mally developing turbulent channel fiow laden with par- 
ticles large enough to neglect Brownian diffusion (which 



becomes important only for particle diameters smaller 
than 1 fim) but small enough to ensure stability against 
the inertia-dominated non-homogeneous distribution [12] 
and consequent near- wall accumulation |10I13) . 



II. METHODOLOGY 



A. Flow field equations 



With reference to the schematics of Fig. [U particles are 
introduced in a turbulent channel flow with heat trans- 
fer. Assuming that the fluid is incompressible and New- 
tonian, the governing balance equations for the fluid (in 
dimensionless form) read as: 
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where Ui is the i component of the velocity vector, p 
is the fluctuating kinematic pressure, Si,i is the mean 
pressure gradient that drives the flow, T is the tempera- 
ture. Re* is the shear (or friction) Reynolds number and 
Pr is the Prandtl number. The shear Reynolds number 
is deflned as Re* = u*h/iy, based on the shear veloc- 
ity, u*, on the half channel height, h, and on the fluid 
kinematic viscosity, v. The shear velocity is defined as 
u* = (tuj/p)^/^, where is the mean shear stress at 
the wall and p is the fluid density. The Prandtl num- 
ber is deflned as Pr = jiCp/k where p, Cp and k are 
the dynamic viscosity, the specific heat and the ther- 
mal conductivity of the fiuid, respectively. All variables 
considered in this study are reported in dimensionless 
form, represented by the superscript +, which has been 
dropped from Eqns. ([l]) to ([3]) for ease of reading, and 
expressed in wall units. Wall units are obtained by tak- 
ing u* , V and the shear (or friction) temperature T* as 
the reference quantities employed for normalization. The 
shear temperature is defined as T* — q^/pcpU* where 
Qw = k- VT!u, is the mean heat flux at the wall and VT„, is 
the wall-normal component of the temperature gradient 
at the wall. In Eqns. ^ and ^ the momentum-coupHng 
term f2w and the energy-coupling term g2io are defined 
in terms of momentum and energy flux per unit mass, 
namely f2w = F2w/mp and q2w = Q2w/{cp ■ rup) where 
TOp is the particle mass. These terms are introduced to 
model, as point sources, the influence of the particles on 
the fluid velocity and temperature flelds (two-way cou- 
pling approach). More details on these terms are given 
in Sec. Ill CI The reference geometry consists of two 
infinite fiat parallel walls; the origin of the coordinate 
system is located at the center of the channel and the 
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Figure 1: Sketch of the computational domain. 



, y~ and z— axes point in the streamwise, spanwise 
and wall-normal directions respectively (see Fig. [J). The 
calculations are performed on a computational domain 
of size Anh x 2'Kh x 2h in x, y and z, respectively. Pe- 
riodic boundary conditions are imposed on both velocity 
and temperature in the homogeneous x and y directions; 
at the wall, no-slip condition is enforced for the momen- 
tum equation whereas constant temperature condition is 
adopted for the energy equation. Specifically, the tem- 
perature on the boundaries is held constant at the uni- 
form values < T'^ >= 80° C for the hot wall (source 
of heat) and < >~ 20° C for the cold wall (sink of 
heat), where <> denotes averaged values. The condition 
on temperature has been chosen because our aim is to 
simulate the problem of thermally-developing forced con- 
vection in a channel where supply of heat from a source 
and release of heat to a sink are considered as constant 
temperature processes. 



B. Particle equations 

Large samples of heavy particles with diameter dp = 
4 /im and 8 /im and with density pp = 19.3 • 10^ kg 
(gold in water) are injected into the flow at concentration 
high enough to have significant two-way coupling effects 
in both the momentum and energy equations but negligi- 
ble particle-particle interactions. The particle dynamics 
is described by a set of ordinary differential equations for 
position, velocity and temperature. For particles heav- 
ier than the fiuid {pp/p ~ 20, as in the present case), 
Armenio and Fiorotto [22] showed that the most signif- 
icant force is Stokes drag. Other forces acting on the 
particles, such as hydrostatic force, Magnus effect. Bas- 
set history force and added mass force are not taken into 
account since they are assumed to be negligible (orders of 
magnitude smaller) because of the specific set of physical 
parameters of our simulations |22|23j . We also neglected 



the Brownian force, which is proportional to 1/dp and 
becomes important for particle diameters less than 1 
[24] . With the above assumptions, a simplified version of 
the Basset-Boussinesq-Oseen equation for particle mo- 
mentum balance [25] is obtained. In vector form: 
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where Xp is the particle position, Up is the particle ve- 
locity, u is the fiuid velocity at particle position and Cd 
is the drag coefiicient. The formulation of the drag coef- 
ficient Cd follows the non-linear approximation reported 
in Schiller and Naumann [26l: 
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where Rcp = dp | Up — u | /j/ is the particle Reynolds 
number based on the relative particle-to-fiuid velocity. 
The correction for Cd is necessary when Rcp does not 
remain small {Rep > 1). 

In this study, the particles exchange both momentum 
and heat with the carrier fiuid. The equation describing 
the temperature evolution of the dispersed phase can be 
derived from the energy balance on a particle, written un- 
der the assumption of convective heat transfer occurring 
through the particle surface (the contribution of radiation 
is thus neglected). The inter-phase heat transfer rate of 
a spherical particle in motion relative to the surrounding 
fiuid can be written as: 



Qc = Nu ■ ndp ■ kf ■ (Tg - T) , 



(7) 



where Nu is the Nusselt number, kf is the thermal 
conductivity of the fiuid, T is the temperature of the 
fiuid and Tg is the temperature at the particle surface. 
The Nusselt number is given by the well-known Ranz- 
Marshall correlation 1271: 



Nu = 2 + 0.6 • Re^^ ■ Pr^^ 



(8) 



and accounts for changes of the heat transfer rate due 
to the relative motion between the particle and the sur- 
rounding fiuid. Rearranging Eqn. one can write: 
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where T/ is the temperature of the fiuid at particle po- 
sition, Tp is the temperature of the particle and tt ~ 
{cpPpdp^)/{12kj:) is the particle thermal response time. 
In this study, we considered uniform temperature inside 
the particle, namely Tp = T^: this assumption is justified 
for particles with Biot number smaller than 0.1 [25]. The 
particle Biot number is defines as: 
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where hp and kp are the convective heat transfer coeffi- 
cient and the thermal conductivity of the particle. In the 
present simulations, the particle Biot number is C(10^^). 

C. Modeling of two-way coupling 

The source terms f2w and q2w in Eqns. ^ and ([3]) 
arise because of the momentum transfer due to the drag 
force on the particle and because of the convective heat 
transfer to/from the particle, respectively. Calculation of 
these coupling terms is done applying the action-reaction 
principle to a generic volume of fluid O containing a par- 
ticle. Focusing on the f2w term (the extension for the 
term q2w is straightforward), we have: 

/ f2»(x)dr! = -f/i„ , (11) 
Jn 

where f/i„ is the force exerted by the fluid on the parti- 
cle. The term £2^,, which represents the feedback of the 
dispersed phase on the fluid, can be obtained by adding 
the contributions of each particle: 

Up 

f2»=E(f2™)' (12) 

where Up is the total number of particles. The evalua- 
tion of each contribution ff^ is obtained using the point 
source approximation |28I29| : 

(L = -{fiuS{^ - xp) , (13) 
where (5(x — Xp) is the Dirac delta function. 

D. DNS methodology 

In this study a DNS of fully-developed channel flow 
with heat transfer is performed. The governing equa- 
tions for the fluid, Eqns. (l)-(3), are discretized using a 
pseudo-spectral method based on transforming the fleld 
variables into wavenumber space, using Fourier represen- 
tations for the periodic (homogeneous) directions and 
a Chebyshev representation for the wall-normal (non- 
homogeneous) direction. As commonly done in pseu- 
dospectral methods, the convective non- linear terms are 
flrst computed in the physical space and then trans- 
formed in the wavenumber space using a de-aliasing 
procedure based on the 2/3- rule; derivatives are eval- 
uated directly in the wavenumber space to maintain 
spectral accuracy. Time advancement of the equations 
is performed using an explicit two-stage Euler/Adams- 
Bashforth scheme for convective terms and an implicit 
Crank-Nicolson method for the viscous terms. The time 
step used is dt~^ = 0.045 in wall units. More details on 
the numerical scheme can be found in [30] . 

DNS calculations were performed at Re* = 150, cor- 
responding to a bulk (average) Reynolds number Rcb = 



Ubh/v = 1900, where Ub is the bulk velocity. At this 
Reynolds number, two different values of the Prandtl 
number were considered (see Table [J). First, a simula- 
tion at Pr = 0.71 (run Rl in Table H]) was performed to 
validate the flow solver against numerical data available 
in the literature [H] for the situation of turbulent flow 
of air (p = 1.3 kg m~^, v = 15.7 • 10~^ w? j s) at ambi- 
ent temperature in a channel of half height h = 0.02 m. 
In this simulation, we have u* = 0.11775 m and 
Uh = 1.49 TO . Second, simulations at Pr = 3 (corre- 
sponding to runs i?2 to i?4 in Table|l| has been performed 
to study heat transfer modiflcations in a particle-laden 
turbulent flow of water {p = 10'^ kg m~^, v = 10~^ m? / s) 
at the temperature < T >= 50° C in a channel of half 
height h = 500 /ito. The temperature < T > corresponds 
to the mean temperature between < > and T <^>. 
In this case, u* = 0.3 to and Ub = 3.8 m s^^. 

The computational domain has dimensions 1885 x 
942 X 300 in wall units and was discretized using an Eu- 
lerian grid made of 128 x 128 x 129 nodes (correspond- 
ing to 128 x 128 Fourier modes and to 129 Chebyshev 
coefficients in the wavenumber space). The grid spac- 
ing is uniform in the homogeneous directions and corre- 
sponds to spatial resolutions equal to Ax"*" = 14.72 and 
Ay+ — 7.35; the nodes along the wall-normal direction 
are clustered near the wall corresponding to spatial res- 
olutions from Az+ = 0.0452 at the wall to Az+ = 3.68 
at the centerHne. The wall-normal grid spacing is al- 
ways smaller than the smallest local flow scale and, thus, 
it fulfllls the requirements imposed by the point-particle 
approach (see Sec. HE). 



E. Lagrangian particle tracking 

The complete set of equations which describes the time 
evolution of particle position, velocity and temperature 
in the turbulent flow fleld is given by Eqns. (|4]), ^ and 
([9]). To solve for these equations, we have coupled the 
DNS flow solver to a Lagrangian tracking routine. The 
routine uses G'^'-order Lagrangian polynomials to inter- 
polate the fluid velocity components and the fluid tem- 
perature at particle position. The performance of the 
interpolation scheme is comparable to that of spectral 
direct summation and to that of an hybrid scheme which 
exploits 6*'*-order Lagrangian polynomials in the stream- 
wise and spanwise directions and Chebyshev summation 



RUN Re* u* Pr dp d+ Tp St tt Str 

[m/s] [jim] [lis] (= r+) [jis] (= r+) 

-Rl 150 0.11775 0.71 — — — — — — 

-R2 150 0.3 3 — — — — — — 

-R3 150 0.3 3 4 1.2 17.16 1.56 5.49 0.5 

i?4 150 0.3 3 8 2.4 68.62 6.24 2.196 2.0 



Table I: Summary of the simulation parameters. 
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in the wall-normal direction. A 4*^-order Runge-Kutta 
scheme is used for time advancement of the particle equa- 
tions. The timestep size is equal to that used for the fluid 
{6t+ = 0.045). The total tracking time was t+ = 2000 
for the 4 fim particles and t'^ = 1500 for the 8 fim par- 
ticles. We remark here that these simulation times are 
not sufficient to achieve a statistically steady state for 
the particle concentration, yet they are long enough to 
obtain converged velocity and temperature statistics and 
to highlight qualitatively the effect of the particles on 
the heat transfer rate. Particles are treated as pointwise, 
rigid spheres (point-particle approach) and are injected 
into the flow at average mass fraction, high enough 
to have a two-way coupling between the particles and the 
fluid ($„i - 10-2) |32I33) . Possible effects due to inter- 
particle collisions are neglected. At the beginning of the 
simulation, particles are distributed randomly over the 
computational domain and their initial velocity and tem- 
perature are set equal to those of the fluid at the particle 
initial position. Periodic boundary conditions are im- 
posed on particles moving outside the computational do- 
main in the homogeneous directions, whereas perfectly- 
elastic collisions at the smooth walls are assumed when 
the particle center is less than one particle radius away 
from the wall. No speciflc boundary condition is needed 
for the particle temperature equation (Eqn. ^ since the 
integration of this equation follows the integration of the 
particle momentum equation (Eqn. [5]) and it requires 
only the knowledge of the initial condition. In the sim- 
ulations presented here, large samples of 800, 000 parti- 
cles have been considered for each value of Re* and Pr. 
We remark here that tracking of O (10^) particles two- 
way coupled with the fluid requires a huge computational 
effort in terms of both computational cost of the simu- 
lation and disk storage availability, considering also the 
rather long tracking times achieved in the simulations. 
Each particle sample is characterized by different val- 
ues of the particle response times. Table [J summarizes 
the complete set of parameters relevant to the simula- 
tions of particle dispersion, including the particle Stokes 
numbers, St and Stx- The particle Stokes number cor- 
responds to the non-dimensional particle response time 
and is obtained using the viscous timescale Tf = i^/u'^ as 
reference. In the present study, we have St = — Tp/rf 
and Str = r^t = Tx/Tf. 

III. RESULTS AND DISCUSSION 

A. Unladen turbulent channel flow with heat 
transfer 

In this paragraph we examine the statistics relative to 
velocity and thermal variables for the base case of un- 
laden fluid. We examine the results relative to the sim- 
ulations Rl and i?2, performed at the same Reynolds 
number {Re* = 150) and at Prandtl numbers equal to 
0.71 and = 3, respectively. Velocity and temperature 



statistics will be examined and compared against litera- 
ture reference cases |19|31j . 



1. Velocity field 

The mean velocity proflle will not be shown since it col- 
lapses onto the logarithmic law of the wall perfectly [13] 
and matches previous results obtained in Refs. |19|31j . 
In Fig. ^ the root mean square (rms) of the fluctuations 
of the velocity components, < ufrmsi^'^) plotted 
as a function of the wall normal coordinate in wall units, 
z+, and compared against the results of Kasagi and lida 
[31] . The agreement is generally good showing small dif- 
ferences which may be due to marginal statistical sam- 
pling of the time series. In the two simulations relative to 
the base fluid only, the velocity fleld depends only on the 
pressure gradient (namely on the shear Reynolds num- 
ber. Re*) and it is not influenced by the value of the 
Prandtl number - forced convection. 



2. Temperature field 

The behavior of the fluid temperature averaged over 
the homogeneous directions {x and y) is shown in Fig. 
([3]) for the two values of the Prandtl number, Pr = 0.71 
and Pr = 3 respectively. In Fig. ((3^), the temperature 
is made dimensionless in outer units, indicated by the 
superscript -, as follows: 



< T- >= 



<T> -T„ 



AT 

z_i J. 



(14) 



where = {Th + Tc)/2 is the average centerline tem- 
perature and ATm = {Th — Tc)/2 is the temperature 
difference between the walls. The average temperature is 
shown as a function of the wall normal coordinate made 
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Figure 2: Rms of fluid velocity components, < 'u-trmsi^^) 
for single-phase (unladen) turbulent channel flow. 
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Figure 3: Mean fluid temperature at Pr = 0.71 and Pr = 3. 
Panels: (a) < T~{z~) > (linear-linear scale), (b) < T^{z'^) > 
(log- linear scale). 



dimensionless by the half channel height, h. The results 
of the present work are compared and assessed against 
those computed by Kasagi and lida [31] for Pr = 0.71 
and against those computed by Na et al. [H] for Pr ~ 3. 
Results match perfectly showing that the simulation pa- 
rameters - i.e. resolution, average time, etc. - are ap- 
propriate for the problem under investigation. 

The same averaged temperature profiles are plotted 
in Fig. ([3b) using a semi-logarithmic scale just to show 
possible minor differences among the current profiles and 
the benchmark profiles |19|31j . Temperature profiles are 
shown as a function of the wall normal distance, both 
being expressed in wall units. In this case, the tem- 
perature is defined as the difference between the local 
average temperature and the wall temperature (both 
relative to the average centerline temperature) normal- 
ized by the friction temperature, T* = qw/{pfCpU*): 
< T+ >= (< T > -Tw)/T*. 

The computed data are also compared with standard 
analytical correlations used to estimate the temperature 
wall dependence. These correlations have the following 



< T+ >= 



< T > -T, 



w 



Pr ■ z+ 

^lnz+ 



Bo 



if z+ < 11.6 
if z+ > 11.6 
(15) 

where ke = 0.21 and Be = -4.4 and 14.7 for Pr = 0.71 
and for Pr = 3, respectively. 

As broadly known, the profiles show the existence of 
a (near-wall) diffusive sublayer, the thickness of which 
varies with the Prandtl number and is approximately 
equal to M+ = 6.5 for Pr = 0.71 and = 4.5 for 

Pr = 3. In Figs. ((SK) and ([3t>) we can further observe 
that the temperature gradient at the wall is strongly de- 
pendent on the value of the Prandtl number. This be- 
havior is shown by Eqn. ifTSj) . from which it is clear that 
the temperature in the viscous sublayer depends linearly 
on the Prandtl number. 

In problems which involve turbulent heat transfer, the 
Prandtl number is also important to estabHsh the small- 
est spatial scale for the temperature field, rjg, which can 
be expressed |34|35j as a function of the Kolmogorov 
scale, ?7fe, as follows: 



1 



Pr) 



3/4 



for Pr < 1, and: 



Pr 



1/2 



(16) 



(17) 



for Pr > 1. The previous relations confirm that for a 
given value of rjk (we remind here that the average value 
of rjk depends only on the Reynolds number, which is 
equal to 150 in the present simulations), the smallest tem- 
perature scale decreases for increasing Prandtl numbers. 

The behavior of rms of the temperature field fiuctu- 
ations, < r+„s(z+) >, is shown in Fig. Results 
obtained by Kasagi and lida [31] for Pr = 0.71 and by 
Na et al. [19] for Pr = 3 are also shown for benchmark 
and comparison. Results are presented in wall variables. 
Focusing on the Pr = 0.71 case, we observe that the 
temperature intensity reaches a maximum at the channel 
centerline, not mimicking the behavior of the fiuctuations 
of the velocity field which all reach their peak in the wall 
proximity (see Fig. This was discussed by Lyons 

et al. [H], who attributed this different behavior to the 
temperature boundary conditions which force a non-zero 
temperature gradient at the center of the channel. In the 
Pr = 3 case, the temperature intensity reaches a maxi- 
mum in the near-wall region. This observation is related 
to the Prandtl number effect on < T^^{z) >, indicating 
that the range of wavenumbers in the thermal fluctuating 
field increases with Pr, for which the spectral functions 
of the velocity fields are negligible. As can be also ob- 
served, the increase in the Prandtl number corresponds 
to a shift of the peak value of the temperature fiuctua- 
tions toward the wall. Both in the current results and in 
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Figure 4: Rms of fluid temperature, < T^rnsi^^) >, at Pr 
0.71 and at Pr = 3. 



those by Kasagi and lida [3T] for Pr = 0.71, the relative 
peak of the temperature fluctuations is located around 
= 20. The location of the peak for the higher Prandtl 
number {Pr = 3) has moved to z"'" = 9, roughly. Small 
discrepancies due perhaps to the marginally statistical 
sample are observed between our results and those of Na 
et al. fT9l. 



B. Particle-laden turbulent channel flow with heat 
transfer and momentum/energy two way coupling 

The main purpose of this section is to analyze the mod- 
ifications on the flow field due to the mutual interaction 
between fluid and particles. In particular, results ob- 
tained from two-way coupling simulations at Pr = 3 are 
compared against those obtained from the correspond- 
ing one-way coupHng simulations, in which particle feed- 
back on the flow is neglected: f2w = and q2w = in 
Eqns. ([2|) and ([3]). We remark here that, even though 
the problem of fluid-particle momentum coupling in two- 
phase flows has been widely investigated |lll28l33j . much 
less effort has been devoted to the energy coupling prob- 
lem. Currently, numerical studies on this problem are 
available for homogeneous shear flow [36] and for homo- 
geneous isotropic turbulence [37]. To our knowledge, this 
is the first attempt to study (by means of DNS) both mo- 
mentum and energy coupling between fluid and particles 
in wall-bounded flow. 



Figure 5: Mean streamwise fluid velocity, < u^{z^) >, 
at Pr = 3: comparison between unladen flow (no momen- 
tum/energy coupling) and flow laden with dp = 4 fj,m parti- 
cles (with momentum/energy coupling). 



streamwise and spanwise directions) and time (over a 
time span of 180 wall units). As expected, mean velocity 
profiles, normaHzed here by the unladen fiow shear veloc- 
ity u* , deviate only slightly (if not negligibly, as in the 
viscous sublayer) from each other. A careful examination 
of Fig. ^ indicates that the effect of particles is to shift 
the velocity profiles slightly toward smaller values in the 
buffer region (5 < z"*" < 30) and toward higher values 
in the outer region (z"*" > 30). Comparison of the mean 
velocity profiles for the 8 /im particle case (not shown 
here for brevity) indicate no observable effect. 

The behavior of the turbulence intensities, given by the 
rms of the fluid velocity fluctuations, < uf^^g{z^) >, 
shows larger differences as presented in Fig. Again, 
the soHd lines refer to the simulation without particles, 
whereas the dot-dashed and the dashed lines refer to the 
two-way coupling simulation with the 4 fim particles and 
with the 8 fim particles, respectively. It appears that par- 
ticles do not affect much the intensities in the near-wall 
region but do substantially change them in the buffer re- 
gion, particularly where the proflles develop a peak, and 
at the channel centerline. For each rms component, lo- 
cal changes of opposite sign depending on particle inertia 
are observed with respect to the reference unladen-flow 
values. 



2. Temperature field modifications by particles 



1. Velocity field modifications by particles 

The effect of particles on the mean fluid velocity, 
< u^{z'^) >, is shown in Fig. |[5|). The solid line refers 
to the simulation without particles, while the dot-dashed 
line refers to the two-way coupling simulation with the 
4 particles. Proflles are averaged in space (over the 



Fig. ^ shows the mean fluid temperature profiles in 
inner units, < T~[z~) >, (Fig. [TK) and in wall units, 
< T+{z+) > (Fig. [7b). Lines are as in Fig. Visual 
inspection of Fig. |[7K) does not reveal significant differ- 
ences in the profiles. However, computing the local value 
of the wall-normal temperature gradient, d<T(z)>/dz, 
right at the wall, we obtain an increase of roughly 2 % 
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Figure 6: Rms of mean fluid velocity components, < 
utrmsi^~^) >i = 3: comparison between unladen 

flow (no momentum/energy coupling) and flow laden with 
dp — 4 iim particles and with dp = 8 /im particles (with 
momentum/energy coupling). Panels: (a) streamwise rms, 
< i'-x,rms{'^^) >; (b) spanwise rms, < iiy,rms(-z^) >; (c) wall- 
normal rms, < uf {z+) >. 



for the two-way coupling simulation with the 4 
particles and a decrease of 0.2 % for the two-way cou- 
pling simulation with the 8 particles. This finding 
is important because changes produced by the particles 
to the wall-normal temperature gradient directly reflect 
upon the heat flux at the wall, q^, through the following 
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Figure 7: Mean fluid temperature at Pr = 3: comparison 
between unladen flow (no momentum/energy coupling) and 
flow laden flow with dp — 4 urn particles and with dp — 8 
/im particles (with momentum/energy coupling). Panels: (a) 
< T-{z-) > (lin-lin scale), (b) < T+{z+) > (log-lin scale). 



expression: 



1 d < T(z) > 
Pr dz 



(18) 



In turn, a change in the value of Qw will eventually cor- 
respond to a change in the value of the total turbulent 
heat flux qz.tot, deflned as: 



qz.tot 



< T'w' > , 



(19) 



under fully developed conditions. Besides being perhaps 
the most signiflcant result of the present paper from a 
quantitative viewpoint, the increase of removable heat 
flux at the wall for the smaller 4 fim particles case and 
the decrease of removable heat flux for the larger 8 fim 
particles case have also an effect on the friction tempera- 
ture, T*. According to its deflnition (see Sec. IIII A|) . the 
friction temperature computed in the two-way coupling 
simulations with the 4 fim particles will increase with re- 
spect to the unladen flow case, whereas it will decrease 
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Figure 8: Rms of fluid temperature, < TXns{z^) >, at Pr — 
3: comparison between unladen flow (no momentum/energy 
coupling) and flow laden with dp = 4 /im particles and with 
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Figure 9: Mean particle streamwise velocity, < u^^p{z^) >, at 
Pr = 3. Symbols: (O) dp — A /im particles, (□) dp = 8 fj,m 
particles. The mean fluid streamwise velocity proflle relative 
to the unladen flow case (solid line) is also shown for sake of 
comparison. 



in the simulations with the 8 fim particles. As a con- 
sequence, starting from the reference profile relative to 
the unladen fiow case, the mean fiuid temperature pro- 
file computed for the 4 particle case will shift toward 
smaller values (particularly outside the viscous sublayer) 
whereas it will shift toward slightly higher values when 
the larger 8 fim particles are considered, as apparent from 

Fig. m- 

Finally, the rms of the fiuid temperature fiuctuations, 

< T+^^{z+) >, are shown in Fig. ([8]). The 4 urn 
particles produce a slight decrease in the peak value of 

< T+„^(2+) > and an increase outside the buffer region; 
smaller modifications (mostly limited to the core region 
of the fiow) are produced by the larger 8 /xto particles. 



3. Influence of particle inertia and of particle thermal 
inertia 

The interactions between particles, turbulent momen- 
tum transport and turbulent heat transport are infiu- 
enced by the particle response times. In this section we 
analyze the Eulerian statistics of the particles in com- 
parison with those of the fiuid. Specifically, particles will 
acquire and lose momentum and heat at a rate propor- 
tional to the inverse of their response times, Tp and tt 
respectively. 

Fig. ^ shows the mean streamwise velocity profile, 
< u^^plz'^) >, averaged in time and along the homoge- 
neous directions, for both fiuid and particles as a function 
of the wall-normal coordinate, z+. Differences are readily 
visible, with a consistent evidence of the effect of particle 
inertia. Smaller particles {St = 1.5, dp = 4 /im) behave 
more like fiuid tracers and their average velocity profile 
(circles) almost match that of the fiuid (solid line). As 
the particle response time increases {St = 6, dp — 8 tim), 



the average particle velocity (squares) is seen to lag the 
average fiuid velocity, in particular outside the viscous 
sublayer. Present results agree well with those reported 
by van Haarlem et al. |38j and by Portela et al. [39], 
who however used a one-way approach. This shows that, 
for the current size, density and overall concentration 
of the particles, modifications to the mean streamwise 
velocity profile due to momentum and energy coupling 
appear negligible. As already observed by van Haarlem 
et al. [38], inertial particles dispersed in a turbulent fiow 
do not sample the fiow field homogeneously and tend to 
avoid areas of high vorticity preferring areas character- 
ized by lower-than-mean streamwise fiuid velocity and by 
high strain rate. This gives the characteristic velocity lag 
in the region 5 < z"*" < 50. This effect is confirmed by 
data shown in Fig. (fTOk ). where the root mean square 
of particles streamwise velocity (circles and squares) is 
compared to that of the fiuid (solid line) . Fluctuations of 
particle streamwise velocity are larger than those of the 
fiuid, this difference becoming more evident as particle 
response time increases [39]. From a physical viewpoint, 
the difference in the streamwise values suggests that the 
gradients in the mean fiuid velocity can produce signifi- 
cant fiuctuations of the streamwise particle velocity. This 
effect seems more pronounced in the case of heavy parti- 
cles, with larger Stokes number and a longer "memory". 
An opposite behavior is observed in the spanwise direc- 
tion and in the wall-normal direction (Figs. (fTOb ) and 
(flOh) , respectively) , where the fluid velocity fleld has zero 
mean gradient. The turbulence intensity of the particles 
is very close to that of the fluid for particles with small 
inertia {St =1.5, circles) and lower than that of the fluid 
for particles with higher inertia {St = 6, squares). This 
is mainly due to two mechanisms acting in tandem. The 
flrst mechanism is preferential concentration of particles 
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Figure 10: Rms of particle velocity components at Pr = 3. 
Symbols: (O) 4 /i?n particles, (□) 8 /^m particles. Panels: 
(a) streamwise rms, < u^^p^^msiz'^) >, (b) spanwise rms, < 
Uy,p,rmsiz'^) >, (c) wall-normal rms, < uJ,p,rms(-2^) >• The 
rms of fluid velocities relative to the unladen flow case (solid 
line) is also shown for sake of comparison. 



in low-speed regions, characterized by lower turbulence 
intensity [13] . The second is the filtering of high frequen- 
cies or wavenumber fiuctuations done by particles due to 
their inertia. The inertial filtering damps turbulence in- 
tensity of the particles in the wall-normal direction and 
in the spanwise direction. Similar filtering effects have 



V 





Figure 11: Mean particle temperature at Pr = 3. Sym- 
bols: (O) 4 /im particles, (□) 8 ^tm particles. Panels: (a) 
< Tp{z^) > (linear-linear scale), (b) < Tp{z~^) > (log-linear 
scale). The mean fluid velocity relative to the unladen flow 
case (solid line) is also shown for sake of comparison. 



been observed in homogeneous turbulence [40] . 

The statistical behavior of the thermal field for the 
dispersed phase is presented in Figs. ifTTj) and lfT2|). Fig. 
(fTTj) shows the mean temperature profiles for particles 
and fiuid as a function of the wall-normal coordinate. In 
Fig. ifTTk ) variables are expressed in outer units, whereas 
in Fig. (fTTb ) variables are shown in wall units and in 
semi-logarithmic scale. The solid line refers to the aver- 
age unladen fiuid temperature, while circles and squares 
refer respectively to the simulations with the 4 fim par- 
ticles and with the 8 fim particles. Fig. ifTTk) indi- 
cates that, for both particle sets, the particle temper- 
ature in the near-wall region (0 < z+ < 30) is higher 
than that of the fiuid, while it reaches lower values in 
the outer region (30 < 2+ < 150). Particle temperature 
higher than that of the fiuid is probably due to the de- 
position/resuspension mechanisms, which bring particles 
to the near-wall region (where they are characterized by 
high temperature value differences) and drive them to- 
ward the core region [10|13j . When a particle initially 
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Figure 12: Rms of particle temperature, < Tj^^rms{z^) >, at 
Pr = 3. Symbols: (O) 4 nm particles, (□) 8 particles. 
The rms of fluid temperature relative to the unladen flow case 
(solid line) is also shown for sake of comparison. 



close to the wall is entrained toward the core region, it 
will adapt its temperature at a rate depending on particle 
thermal inertia, thus maintaining a temperature higher 
than that of the fluid. In the outer region, the perspective 
is overturned, and particles are characterized by a tem- 
perature lower than that of the fluid. The effect of par- 
ticle thermal inertia is also visible by observing the rms 
of particle temperature proflles, < >, shown 

in Fig. (fT2| . These proflles are qualitatively similar to 
those for the streamwise rms component (Fig. \Wk). due 
to the role played by particle thermal inertia in trans- 
ferring heat, similar to that played by particle inertia in 
transferring momentum. Although a qualitative similar- 
ity between the temperature and the streamwise velocity 
rms can be observed, quantitative differences are visible. 
In particular, the peak of the streamwise velocity rms 
exhibits values higher than the temperature rms. These 
quantitative differences are due to the fact that particle 
inertia is larger than particle thermal inertia and causes 
streamwise velocity fluctuations higher than thermal fluc- 
tuations. 



4- Instantaneous features of turbulent transfer of 
momentum, heat and particles 

A number of DNS-based works have been successful in 
clarifying the role of the instantaneous realizations of the 
Reynolds stresses in transferring momentum [H], heat 
[19] and particles [lO]. We refer to these works and to 
the references therein for a more comprehensive descrip- 
tion of the phenomena. It suffices here to show on a 
qualitative basis some relations among momentum, heat 
and particles transfer mechanisms. 

To this aim, we complement the statistical analysis of 
the previous sections by showing the complex interac- 
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Figure 13: Front view of particle instantaneous distribution 
superimposed to the fluid velocity field (top) and to the fiuid 
temperature field (bottom) at Pr = 3. 



tions among velocity field, temperature field and parti- 
cles distribution. In particular, we examine simulation 
i?3 {Re* — 150, Pr = 3 and dp = 4 fim particles). We 
consider a cross section of the flow fleld, perpendicular 
to the mean velocity, cut along a, x — y plane. Consider 
flrst Fig. ifTSk ): color isocontours correspond to values 
of the streamwise fluid velocity whereas circles represent 
particles (drawn larger than the real size for visualization 
purposes) colored by their wall-normal velocity. In this 
Figure, green particles have wall-normal velocity directed 
toward the walls, while black particles have wall-normal 
velocity directed away from the walls. In this way, just 
by analyzing particle wall-normal velocity, it is possible 
to detect particle transfer fluxes toward the wall and par- 
ticle transfer fluxes away from the wall. It is possible to 
observe that fluxes of particles are associated with fluxes 
of streamwise momentum (called sweeps if directed to- 
ward the wall and ejections if directed away from the 
wall, as discussed in several previous papers |10|13) ). 

To draw a link between the heat transfer mechanisms 
and the momentum transfer mechanisms we consider 
also Fig. ifTSb). where the particles are shown super- 
imposed to the temperature fleld for the same section of 
Fig. (fT3k ). In this case color isocontours are the values 
of the fluid temperature, whereas circles represent par- 
ticles colored by their temperature: speciflcally, green 
particles have higher-than-mean temperature, blue par- 
ticles have lower-than-mean temperature, and red parti- 
cles have temperature close to the mean. The Reynolds 
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Figure 14: Front view of particle instantaneous distribution 
(left panels) and corresponding particle wall-normal concen- 
tration, C/Co (right panels) at Pr = 3. Panels: (a-b) 8 fim 
particles, (c-d) 4 jim particles. 



transport analogy is clearly demonstrated by the behav- 
ior of the instantaneous temperature field: focusing on 
the wall at z+ = 0, characterized by higher temperature, 
we observe hot fluid plumes raising in correspondence 
with the raising of low momentum fluid within an ejection 
in Fig. (fT3k ). This just demonstrates the effectiveness 
of the Reynolds stresses in transporting the fluid close to 
the wall towards the center of the channel and viceversa. 
Similar observations can be made for the cold wall at 
z'^ = 300. Considering now, again in Fig. (fT3b ). the be- 
havior of the particles in correlation with the temperature 
fleld, we observe that particles with higher temperature 
are ejected from the lower wall and are directed towards 
the upper wall. The opposite occurs to the colder parti- 
cles. In this case, however, it is important to remind that 
particle temperature is not strictly correlated with fluid 
temperature. This behavior can be explained by consid- 
ering that particle trajectories depend on fluid velocity 
and particle inertia, but not on fluid or particle temper- 
ature. Particles, driven by the fluid vortices, can thus 
reach regions characterized by fluid temperature quite 
different from particle temperature, causing an apprecia- 
ble heat exchange between the two phases. 

The efficiency of the overall heat exchange process is 
conditioned by the degree of non-homogeneity of particle 
distribution and by the rate of particle accumulation at 
the wall. This observation is qualitatively corroborated 
by Fig. lfT4|) , where the instantaneous cross-sectional dis- 
tribution (left-hand side panels) and the wall-normal con- 
centration, C/Cq (right-hand side panels) of both parti- 
cle sets are compared at the same time instant. Note 
that the 4 /im particle distribution, shown in Fig. lfT4h ). 
is the same of Fig. ifTSj) . the particle color code being 
that of Fig. ifTSb ): particle concentration is computed 
as particle number density distribution per unit volume 
normalized by its initial value [10]. It is evident that the 
smaller 4 particles, which produce an increase of heat 
transfer at the wall (as discussed in Sec. IIII B[) . exhibit 
a more persistent stability against non-homogeneous dis- 
tribution and near-wall concentration with respect to the 
larger 8 particles. The larger particles tend to form 
clusters in the core region and accumulate at the wall at 
higher rates acting as an additional thermal resistance 
both between the walls and the fiuid and between the 
fiuid and the particles: as a result of this behavior, the 
heat transfer is reduced. 



IV. CONCLUSIONS AND FUTURE 
DEVELOPMENT 

Heat transfer enhancement is a fascinating subject 
with extremely interesting possibilities for application. 
One option to increase heat transfer is to devise a new 
concept of heat transfer media constituted by a base fiuid 
in which suitably-chosen heat transfer agents, precisely 
micro and nano particles, are injected. In this way, the 
fiuid can be a standard fiuid characterized by simplicity of 



use and well-known properties, like water, and the heat 
transfer agents can be heavy-metal, high-heat-capacity, 
dispersed particles (e.g. copper, gold or platinum). 

Current literature trends show the potentials of such 
heat transfer media, called nanofiuids; yet the compli- 
cacy and the cost of experimental methods make it hard 
to understand the intricacy of the mechanisms which gov- 
ern the dynamics of the turbulent heat transfer among 
the fiuid and the particles. Current open questions range 
from the optimal size of the particles, to the optimal 
concentration and practical solutions in real applications. 
Further complicating effects are represented by the par- 
ticle inertia and the particle thermal inertia, which are 
additional parameters. 

The present study represents a first effort made in the 
frame a broader project on the numerical simulation of 
heat transfer in nanofiuids. Our strategic object is to 
investigate the heat transfer mechanisms in nanofiuids 
and to devise a suitable numerical methodology to anal- 
yse their behavior. In this paper, the DNS of fiuid and 
thermal fields is assessed against literature data (Kasagi 
and lida [31] and Na et al. [H]) for one value of the 
shear Reynolds number. Re* = 150, and two values of the 
Prandtl number, Pr = 0.71 and Pr = 3. Further, prelim- 
inary results obtained from DNS of two-phase solid-liquid 
turbulent fiow in a channel with heat transfer are pre- 
sented. Hydrodynamically fully-developed, thermally- 
developing fiow conditions have been considered to inves- 
tigate on heat transfer modulation produced by the dis- 
persion of micrometer sized particles. Two different sets 
of particles are considered, which are characterized by 
dimensionless inertia response times equal to St = 1.56 



and St = 6.24 and by dimensionless thermal response 
times equal to Stx = 0.5 and Stx = 2, respectively. The 
instantaneous features of the velocity field, the temper- 
ature field and particles dispersion are discussed from a 
qualitative mechanistic viewpoint to analyse the role of 
particles as heat transfer enhancement agents. 

Building on the results presented in this paper, an ex- 
tensive numerical work combined with modeling efforts 
will be carried out to achieve a deeper understanding of 
the reason why nanofluids conduct heat so effectively. For 
instance, the increased surface interaction between the 
fluid and the solid particles at the nanoscale as possible 
explanation to the increased heat transferability will be 
investigated. Indeed, for a given volume of material there 
is a greater number of particles as their size decreases; 
perhaps there is more opportunity for the nanoparti- 
cles to conduct the heat. Furthermore, understanding of 
how the molecules of a base fluid keep nanoparticles sus- 
pended, since nanoparticles are still dramatically larger 
than individual molecules needs to be investigated. In 
this context, the effect of Brownian forces on the kine- 
matics of the nanoparticles should be investigated. It 
would be also of interest to investigate the magnitude of 
the van der Waals forces between the particles and their 
effect on the nanofluids dynamics. These forces are usu- 
ally small, but they become strong (and attractive) when 
the distance between particles becomes of the order of 
tenths of nanometers. 
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